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Abstract 

The concept of tensor rank, introduced in the twenties, has been popularized at the 
beginning of the seventies. This has allowed to carry out Factor Analysis on arrays with 
more than two indices. The generic rank may be seen as an upper bound to the number of 
factors that can be extracted from a given tensor. We explain in this short paper how to 
obtain numerically the generic rank of tensors of arbitrary dimensions, and compare it with 
the rare algebraic results already known at order three. In particular, we examine the cases 
of symmetric tensors, tensors with symmetric matrix slices, or tensors with free entries. 

Resume 

La notion de rang tensoriel, proposee dans les annees vingt, a ete popularisee au debut des 
annees soixante-dix. Ceci a permis de mettre en oeuvre I'Analyse de Facteurs sur des tableaux 
de donnees comportant plus de deux indices. Le rang generique pent etre vu comme une 
borne superieure sur le nombre de facteurs pouvant etre extraits d'un tenseur donne. Nous 
expliquons dans ce court article comment trouver numeriquement le rang generique d'un 
tenseur de dimensions arbitr aires, et le comparons aux quelques rares resultats algebriques 
deja connus a I'ordre trois. Nous examinons notamment les cas des tenseurs symetriques, 
des tenseurs a tranches matricielles symetriques, ou des tenseurs a elements libres. 
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1 Introduction 



Generic ranks, defined in the complex field, have been studied for several decades [TO] [H]. How- 
ever, the value of the generic rank for arbitrary dimensions is not yet known in the unsymmetric 
case, and has been known in the symmetric case only recently [1]. The existence itself of the generic 
rank is not ensured in the real case, and there exist in general several typical ranks (see section 12.11 
for definitions). The typical tensorial rank of three-way arrays over algebraically closed fields has 
been much studied in the context of computational complexity theory. Biirgisser, Clausen and 
Shokrollahi ^ Ch.20] give an overview of general results for various classes of arrays; these results 
have been extended in [1]. The study of tensorial rank over the real field has lagged behind. In 
this paper, generic and typical ranks are discussed for various tensor structures. 

The typical rank of three-way arrays over the real field has been relevant for psychological data 
analysis since Carroll and Chang [3] and Harshman [9| independently proposed a method which 
they christened Candecomp and Parafac, respectively. This Canonical Decomposition (CanD) 
method generalizes Principal Component Analysis to three-way data, by seeking the best least 
squares approximation of the data array by the sum of a limited number of rank-one arrays. In 
2-way analysis, the rank of the data matrix is the maximum number of components that Principal 
Component Analysis can extract. This property generalizes smoothly to three-way data. That is, 
the rank of a three-way array is the maximum number of components that CanD can extract. 
Thus, the study of typical rank of three-way arrays is of great theoretical importance for CanD. 

Although CanD was developed in a psychometric environment, its main area of applications 
has been Chemometrics, e.g. [13]. In addition to straightforward application of CanD, chemo- 
metricians also use TuckerS component analysis [21j quite often. This is a generalized version 
of CanD which decomposes a data array as weighted sum of rank-one arrays, the weights being 
collected in a so-called core array. Typically, the underlying chemometric model dictates that a 
vast majority of specified core elements vanish. Because there exist admissible transformations 
which generate a vast majority of zero elements in arbitrary arrays, we need tools to tell models 
from tautologies. This is were the concept of typical rank has found another realm of application. 
For instance. Ten Berge and Smilde [i9\ have argued that a sparse core hypothesized by Gurden 
et alterae [7| is indeed a model and not a tautology. Their hypothetical core was a 5 x 5 x 3 array 
with only 5 nonzero entries, hence of rank 5 at most. Because 5x5x3 arrays have a typical rank 
of at least 7, it is clear at once that transformations which yield as few as five non-zero elements, 
starting from any randomly generated 5x5x3 array, do not exist. In this way, the typical rank of 
three-way arrays find applications in distinguishing constrained Tucker3 models from tautologies. 

The paper is organized as follows. In section [21 definitions and historical remarks are provided. 
Next, a numerical algorithm is described in section [3l which is able to compute the generic rank 
of any tensor, symmetric or not. Numerical values are reported in section HJ and compared 
to the already known rank values previously obtained by means of algebraic calculations. The 
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consistency of the results confirm tlie validity of the approach, which can yield generic ranks for 
more complicated structures, such as tensors with symmetric matrix slices (sometimes referred to 
as the INDSCAL model), among others. 

2 Generic and Typical Ranks 
2.1 Definitions 

Let T be a L-way array of dimensions A^^, 1 < i < L, with values in a ring TZ. This array always 
admits a decomposition into a sum of outer products as: 

R 

T = J2 u^'^ ® 1^1'^ ® . . . ® 4"^^ (1) 

where itr is a A£ x 1 array, Vr, and ® denotes the tensor product. 

Now consider an array T with values in a field K. Arrays ui^^ may be considered as vec- 
tors of the linear space K^'. Thus, as a combination of tensor products of vectors, T may be 
considered as a tensor. Under a linear change of coordinate system in each space K^*^, defined 
by a matrix A*-^-*, the tensor is represented by another array, obtained by the multi-linear trans- 
form {A*-^'', A*-^-*, . . .A^^^}. Since it is legitimate once a basis has been defined in the space, no 
distinction will be made in the remainder between the tensor and its array representation. 

The rank of a given tensor T (and by extension, of the array defining its coordinates in a 
given basis) is the minimal integer R such that the decomposition ([T]) is exactly satisfied. This 
decomposition is referred to as the tensor Canonical Decomposition (CanD). 

A property is called typical if it holds true on a set of nonzero volume [1] . This supposes that 
some topology has been defined on k^i^^^x.-.a^l. ^j^jg Zariski topology for instance, or 

an Euclidian topology. A property is said to be generic if it is true almost everywhere. In other 
words, a generic property is typical, but the converse is not true. 

Let Ni be L given positive integers. Then the rank of tensors of size Ni x N2 x ■■■ x Ni is 
bounded, and one can make a partition of the tensor space, according to the rank values. One can 
define typical ranks as the ranks that are associated with subsets of nonzero volume in the latter 
partition. If there is a single typical rank, then it may be called the generic rank. 

For instance, there is a single generic rank if the underlying field K is algebraically closed (as 
the field of complex numbers, C) [H] [1]. But there may be several typical ranks if K is the real 
field, M. 
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2.2 Historical remarks 



The study of typical rank over the real field was initiated by Kruskal [TT] [12], who noted that 
2x2x2 arrays had both rank 2 and rank 3 with positive probability. Kruskal also added a 
few typical ranks for small arrays. Franc ^ discussed some more results, including bounds on 
typical rank. Ten Berge and Kiers [T7] gave a first result of some generality, in solving the typical 
rank issue for all two-slices arrays (that is, arrays of format 2 x J x K). These results were 
further generalized in [15], to include all cases where, for I > J > K, I > JK — J. Additional 
miscellaneous results can be found in [1] [15] [16] [20] . 

When Carroll and Chang developed Candecomp, the main applications they had in mind (a 
scalar product fitting problem related to INDSCAL) involves three-way arrays with slices that are 
symmetric in two of the three modes. Ten Berge, Sidiropoulos and Rocci \l8l noted that this form 
of symmetry would affect typical ranks, and examined a number of cases. Quite often, indeed, 
symmetry of slices appears to entail lower typical rank values. On the other hand, there are also 
cases where symmetry of the slices does not affect the typical rank. A partial explanation for 
this can be found in [20]. Ten Berge et al. also noted that symmetric slices are often double 
centered [3], p. 286], which will further reduce the typical rank. That is, when an array has I 
double centered slices of order J x J, it can be reduced toa/x (J— 1) x (J — 1) array, and its 
typical rank will therefore be the same as that of noncentered symmetric / x (J — 1) x (J — 1) 
arrays. A parallel reasoning can be given for nonsymmetrical double centered slices. A rationale 
for double centering slices in the PARAFAC context can be found in [8l p. 239]. It is easy to show 
that the I x J x K array with I double centered slices has the same typical rank as the uncentered 
/ X (J — 1) X (i^ — 1) array. 



3 Computation of Generic Ranks 

The algorithm proposed is directly inspired from |0. Equation ([T]) can be seen as a parametrization 
of tensor T. In fact, given a set of vectors {ui^^ G K^*, 1 < i < L, 1 < r}, consider the mapping 
(f defined from a known subspace Tr of (K^^ x ^ . . . x K^^)^ onto K^i^^.-A^l g^g. 

R 

{ul^^ eTn, 1 <i < L,l <r < R} ® ® • • • ® ^i^^ 

r=l 

Denote = ^{Tr) the image of this mapping. Then the dimension D of its closure Zr is given 
by the rank of the Jacobian of expressed in any fixed basis of "^^^^'^-^l _ j^^ie Jacobian is of 
maximal rank, that is, if its rank equals the number M of free parameters in 7r, then it means 
that Ris a. typical rank. Actually, R will be either the smallest typical rank, or the generic rank. 
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Note that it is always possible to reach the maximal Jacobian rank by increasing the number of 
terms i?, so that the smallest typical rank is always found. 
This result yields the following numerical algorithm. 

• Express formally the parametrized rank-one tensor term in a canonical basis 

• Express formally the gradient of the latter in this basis 

• Draw randomly the parameters according to an absolutely continuous distribution, and 
initialize matrix J with the numerical value of the gradient 

• While rank( J) strictly increases, do: 

— Draw randomly the parameters according to an absolutely continuous distribution, and 
append this new numerical value of the gradient as a new row block in J 

— Compute the new value oi D — rank (J) 

• Compute the dimension of the fiber of solutions as F = M — D, the difference between the 
number of parameters and the dimension of the image Zj^. 

In order to clarify the description of this algorithm, we give now the exact expressions of the 
Jacobian in various cases. 

3.1 Jacobian for 3rd order asymmetric tensors with free entries 

The mapping takes the form below 

r 

{a{e), 6(£), c(£)} T = ^ a{i) ® b{i) ® c{i) 

taking into account the presence of redundancies, the number of parameters in this parametrization 
is M = -R(A^i + A^2 + -^3 — 2). In a canonical basis, T has the coordinate vector: 

r 

e=i 

where we may decide that a, 6, and c are row arrays of dimension A^i, A^2, and N^, respectively, 
and (8) denotes the Kronecker product. Hence, after r iterations, the Jacobian of (p is the r{Ni + 
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N2 + N3) X N1N2N3 matrix: 





® 6^(1) 




c-(l) 


a{lY 






c-(l) 


aiiy 


® b{lY 




Ins 


I-rii 


® ... 

® i"„2 
® ... 




Ins 


I ni 


® b^{r) 






a{rY 


(g) i"„2 






a{rY 


(g) h{rY 




Ins 



The values of the generic rank obtained with this algorithm, called rangj3(Nl,N2,N3), or 
rangj (N,L) for tensors of arbitrary order L and equal dimensions, are reported in tables [H [21 and 

El 



3.2 Jacobian for 3rd order asymmetric tensors with symmetric matrix 
slices 

In this section, consider tensors of size J x J x K, having symmetric J x J matrix slices. Our 
code rgindscal3(J,K) implements the computation of the rank of the Jacobian below, when its 
size increases according to the algorithm described in section [31 

I J ® 6(1) ® c(l) + b{l) ®Ij® c(l) 
6(1) ® 6(1) 

I J ® 6(r) (g) c + 6(r) ®Ij® c(r) 
6(r) (g 6(r) (g Ik 

After r iterations, this matrix is of size r(J + K) x J'^K. The number of parameters in this 
parametrization is M = R{J + K — 1). Values of the generic rank are reported in table [H 



3.3 Jacobian for 3rd order double centered asymmetric tensors 

Now, take again J x J x K tensors with symmetric J x J matrix slices, but assume in addition 
that every row and column in the latter matrix slices are zero-mean. In order to achieve this, it 
is sufficient to generate vectors 6(r) with zero mean; in otehr words, only J — 1 random numbers 
need to be drawn, the last entry of each vector 6(r) being obtained via hj = —J2j=i^j- The 
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Jacobian takes then the expression below: 



■ j-i, 



^ j-i, 



b(l) ® c(l) + 6(1) ® -1] ® c(l) 
6(1) ® 6(1) 

6(r) (g) c(r) + 6(r) (g) -1] (g) c(r) 
6(r) (g) 6(r) (g) /i^- 



where 1 denotes a column of ones of size J — 1. At the rth iteration, this matrix is of size 
r(J + i^' — 1) X J^K. The number of parameters in this parametrization is M = R{J + K — 2). 
Table O reports some numerical values obtained with the code rgindscal2z. 



3.4 Jacobian for symmetric tensors 

In the case of symmetric tensors of dimension and order L, the mapping is defined from K 

L 



NR 



to the space of symmetric tensors 0, or equivalently to with p = {^^^ as: 



where (g) stands for the tensor (outer) product; once a basis is chosen, the tensor product may 
be replaced by a Kronecker product, yielding exactly the same expression. In the case of order-3 
tensors (L = 3) and after r iterations, the Jacobian of takes the following form, somewhat 
simpler than the previous cases: 

In ® a{l) (g) a(l) + a(l) ®In® a{l) + a(l) (g) a(l) ® In 
In (g) a{r) ® a{r) + a{r) ® In ® a{r) + a(r) (g) a{r) ® In 

This matrix is of size rA^xA^^, but we know that its rank cannot exceed {^^'^) = N{N+l){N+2)/Q. 
The number of parameters in this parametrization is M = RN. Numerical values of the generic 
rank obtained with rangjs(N,L) are reported in table [6l 



4 Numerical results 

The available results on unconstrained, slicewise symmetric, and double centered arrays can be 
compared with the numerical values delivered by the computer codes. 

Tensors with free entries. Table [1] reports typical ranks for 2-slice, 3-slice, and 4-slice 
arrays. All known typical rank values are reported in plain face, and coincide with the results 
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Table 1: Typical ranks for 2-slice, 3-slice, and 4-slice unconstrained arrays. Values reported in 
bold correspond to generic ranks computed numerically. Values separated by commas are known 
typical ranks. 







K=2 






K=3 




K= 


=4 




J=2 


J=3 


J=4 


J=3 


J=4 


J=5 


J=4 


J=5 


1=2 


2,3 


3 


4 


3,4 


4 


5 


4,5 


5 


1=3 


3 


3,4 


4 


5 


5 


5,6 


6 


6 


1=4 


4 


4 


4,5 


5,? 


6 


6 


7 


8 


1=5 


4 


5 


5 


5,6 


6 


8 


8 


9 


1=6 


4 


6 


6 


6 


7 


8 


8 


10 


1=7 


4 


6 


7 


7 


7 


9 


9 


10 


1=8 


4 


6 


8 


8 


8,9 


9 


10 


11 


1=9 


4 


6 


8 


9 


9 


9 


10 


12 


1=10 


4 


6 


8 


9 


10 


10 


10 


12 


1=11 


4 


6 


8 


9 


11 


11 


11 


13 


1=12 


4 


6 


8 


9 


12 


12 


12,13 


13 



Table 2: Generic rank of unconstrained arrays of dimension K x K x K . 



K 


2 


3 


4 


5 


6 


7 


8 9 


R 


2 


5 


7 


10 


14 


19 


24 30 



from rang j 3. As can be verified, the smallest of the known typical rank values within a cell 
coincides throughout with the results from rang j 3. For the unknown entries, we insert the results 
from rang j 3 in bold face. These bold face values represent the smallest typical rank. 

We report values of the generic rank of 3-way arrays with equal dimensions in table [2j Kruskal 
[121 p. 9] refers to a "much studied 9x9x9 array whose rank has been bounded between 18 and 
23 but is still unknown". One can observe that the code rangj3(K,K,K) yields a generic rank of 
^ = 30 for i^" = 9, which shows that the array Kruskal refers to is sub-generic. The values shown 
in table [2] can also be compared to those obtained in the symmetric case (see table [6]). 

Now the algorithm can be run on tensors of order higher than 3. To make the presentation of 
the results readable, table [3] reports values of the generic rank obtained for asymmetric tensors with 
equal dimensions, N, and order L, with an algorithm referred to as rangj (N,L). We also indicate 
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Table 3: (top) Generic rank of unconstrained arrays of equal dimensions, A^, and order L. (bottom) 
Number of remaining degrees of freedom; essential uniqueness of the CanD occurs when F = 0. 



R 


N=2 


N=3 


N=4 


N=5 


N=6 


N=7 


N=8 


L = 3 


2 


5 


7 


10 


14 


19 


24 


L = 4 


4 


9 


20 


37 


62 






F 


N=2 


N=3 


N=4 


N=5 


N=6 


N=7 


N=8 


L = 3 





8 


6 


5 


8 


18 


16 


L = 4 


4 





4 


4 


6 







the dimensionality of the fiber of solutions. This number is simply defined as the difference: 

F(iV, L) = R{N, L) {LN - L + 1) - iV^ 

Only those values of dimension and order for which F = have an essentially unique CanD, that 
is, unique up to usual scale and permutation ambiguities. 

Tensors with symmetric matrix slices. Having verified that rangjS and rangj work 
correctly throughout the cases were the generic/typical rank are known, we next turn to the 
I X J X J arrays with I symmetric slices (Table H]). Again, known values coincide with numerical 
ones delivered by the code rgindscalS. We inserted results obtained from rgindscalS alone in 
bold face. As far as can be determined, all results are again in agreement with previously known 
values. 

Tensors with double centered symmetric slices. When the symmetric slices are also 
row-wise (or column-wise, which is the same thing) zero-mean, the code rgindscal2z yielded the 
values reported in table O Note that the generic rank computed by rgindscal2z(J,I) is the 
same as that computed by tgindscal3(J-l , I), at least according to the values explored in table 
m This confirms the conjecture we made earlier in this paper. 

Symmetric tensors. Finally, we also report in table [H] values obtained with 3- way or 4-way 
symmetric arrays, obtained with the code rangj s. Note that these results have been already 
reported in [5] . The dimensionality of the fiber of solutions is given by: 

It is interesting to compare the ranks with those of the unsymmetric case, obviously larger, re- 
ported in table [21 In particular, by inspection of the values of F, one can observe that uniqueness 
is again rarely met with generic arrays, but less rarely than in the non-symmeric case. 
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Table 4: Typical ranks for I x J x J arrays, with J x J symmetric slices. Values reported in 
bold correspond to generic ranks computed numerically. Values separated by commas are known 
typical ranks. 
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6 


7 


9 


1=7 


3 
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10 


1=8 
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10 


1=9 
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9,10 


11 


i=10 
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6 


10 


11 



Table 5: Typical ranks for IxJxJ arrays, with J x J symmetric slices having zero-mean columns. 
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J=5 
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1=3 
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1=4 
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3 
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1=5 
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1=7 
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6 
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1=8 
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6 
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i=10 
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3 


6 


10 



10 



Table 6: (top) Generic ranks of symmetric arrays of dimension and order L. (bottom) Number 
of remaining degrees of freedom; essential uniqueness of the CanD occurs when F = 0. 
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N=2 


N=3 


N=4 


N=5 


N=6 


N=7 


N=8 


L = 3 
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4 


5 
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10 


12 


15 


L = 4 


3 


6 


10 


15 


21 


30 


42 


F 


N=2 


N=3 


N=4 


N=5 


N=6 


N=7 


N=8 


L = 3 





2 





5 


4 








L = 4 


1 


3 


5 


5 
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Appendix: Computer codes 

The CO de£l are provided in this section in the form of simple Scilab functions (Scilab is a softaware 
made available by Inria, freely to academic laboratories; its syntax ressembles that of Matlab). 

function [Rbar,D,F]=rangj3(Nl,N2,N3) 
// [Rbar,D,F]=rangj3(Nl,N2,N3) 
// P.Comon, 22 May 2006. 

// Rbar is the generic rank of tensors of order 3 and dimensions NlxN2xN3, 
// and D is the dimension of the whole space of square tensors 
// F is the dimension of the variety of solutions 

// if a norm constraint is imposed on loading vectors in each term, 
// F=F=Rbar*(Nl+N2+N3-2)-D 

// Call: getf ("rangj3.sci") ;stacksize(2000000) ; 
// [Rbar,D,F]=rangj3(Nl,N2,N3) 

// Example: K=4;for 1=2:9, for J=3:5,M(I-1, J-2)=rangj3(I, J,K) ;end;end;M 
Il=eye(Nl,Nl) ;I2=eye(N2,N2) ;I3=eye(N3,N3) ; JAC=[] ; 

rankold=0 ; ranknew=l ; 

while ranknew>rankold, 
a=rand(l,Nl, "normal") ;b=rand(l,N2, "normal") ;c=rand(l,N3, "normal") ; 
JAC= [JAC;kron(Il ,kron(b,c) ) ;kron(a,kron(I2,c) ) ;kron(a,kron(b, 13) )] ; 

temp=ranknew; ranknew=rank( JAC) ; 

printf ( ' ranknew=°/oi\n ' , ranknew) 

rankold=temp ; 

end; 

D=ranknew ; 

^These codes can be downloaded from the URL www. i3s .unice . f r/^pcomon, both in Scilab and Matlab 
formats. 
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[lign , col] =size ( JAC) ; Rbar=lign/ (N1+N2+N3) -1 ; F=Rbar* (N1+N2+N3-2) -D ; 

function [Rbar,D,F]=rangj (N,K) 
// [Rbar,D]=rangj (N,K) 
// P.Comon, 21 May 2006. 

// Rbar is the generic rank of square tensors of order K and dimension N, 
// and D is the dimension of the whole space of square tensors 
// F is the dimension of the variety of solutions 

// if a norm constraint is imposed on K-1 loading vectors in each term, 
// F=Rbar*(K*N-K+l)-D 

// Call: getf ("rangj .sci") ;stacksize(2000000) ; 

// for n=2:5,for k=3:4, [R,D]=rangj (n,k) ;G(n-l,k-2)=R;F(n-l ,k-2)=R*(k*n-k+l)-D;printf ('\n') ;end;enc 

// NB: this implementation limited to K=3 or K=4, and to N<6 for reasons of 

// numerical complexity and memory requirements. 

getf ("/users/comon/Documents/Scilab/fonctions/gsm. sci") 

I=eye(N,N) ; JAC=[] ; 

if K==3, 

rankold=0 ; ranknew=l ; 

while ranknew>rankold, 
a=rand(l ,N, "normal") ;b=rand(l ,N, "normal" ) ; c=rand(l ,N, "normal") ; 
JAC= [JAC ; kr on ( I , kr on (b , c ) ) ; kr on (a , kron ( I , c) ) ; kron (a , kron (b , I ) ) ] ; 

temp=ranknew;ranknew=rank(JAC) ; 

/ / temp=ranknew ; [Q , U] =gsm ( JAC) ; ranknew=rank (U) ; 

printf ( 'ranknew=%i\n' ,ranknew) 

rankold=temp; 

end; 

D=ranknew ; 
elseif K==4, 
rankold=0 ; ranknew=l ; 
while ranknew>rankold, 

a=rand(l ,N, "normal") ;b=rand(l ,N, "normal") ; c=rand(l ,N, "normal") ;d=rand(l ,N, "normal") ; 

JAC= [JAC ; kr on ( I , kr on (b , kr on ( c , d) ) ) ; kr on ( a , kr on ( I , kr on ( c , d) ) ) ; kr on ( a , kron (b , kron ( I , d) ) ) ; kr on (a , ki 
temp=ranknew;ranknew=rank(JAC) ; 
// temp=ranknew; [Q ,U] =gsm( JAC) ; ranknew=rank (U) ; 
printf ( 'ranknew=%i\n' ,ranknew) 
rankold=temp; 
end; 

D=ranknew ; 
else 

Rbar=0;D=0; 
end; 
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[lign, col] =size ( JAC) ; Rbar=lign/N/K-1 ; F=Rbar* (K*N-K+1) -D ; 



function [Rbar , D , F] =rgindscal3 (N2 , N3) 
// [Rbar,D,F]=rgindscal3(N2,N3) 
// P.Comon, 29 June 2006. 

// Rbar is the generic rank of tensors of order 3 and dimensions (N2,N2,N3) 
// with symmetric slices across the two first modes 
// and D is the dimension of the whole tensor space 

// F is the dimension of the variety of solutions 

// if a norm constraint is imposed on K-1 loading vectors in each term, 
// F=Rbar*(N2+N3-l)-D 

// rmin is a lower bound on the generic rank, if known. The default is 1. 

// Call: getf ("rgindscal3.sci") ;stacksize(2000000) ; 

// for 1=2:10, for J=2:5, M(I-1, J-l)=rgindscal3(J,I) ;end;end;M 

I2=eye (N2 , N2) ; I3=eye (N3 , N3) ; JAC= [] ; 

raiikold=0 ; raiiknew=l ; 

while ranknew>rankold, 
b=rand ( 1 , N2) ; c=rand (1 , N3) ; 

JAC= [JAC ; kron (12 , kron (b , c) ) +kron(b , kron (12 , c) ) ; kron (b , kron (b , 13) ) ] ; 
temp=ranknew;ranknew=rank(JAC) ; 
printf ( 'ranknew=%i\n' ,ranknew) 
rankold=temp ; 
end; 

D=ranknew ; 

[lign , col] =size (JAC) ; Rbar=lign/ (N2+N3) -1 ; F=Rbar* (N2+N3-1) -D ; 

function [Rbar , D , F] =rgindscal2z (N , K) 
// [Rbar,D,F]=rgindscal2z(N,K) 
// P.Comon, 1st July 2006. 

// Rbar is the generic rank of tensors of order 3 and dimensions (N,N,K) 
// with sjmmetric slices across the two first modes, and doubly centered 
// and D is the dimension of the whole tensor space 
// F is the dimension of the variety of solutions 

// if a norm constraint is imposed on loading vectors in each term, 
// F=Rbar*(N+K-2)-D 

// Call: getf ("rgindscal2z.sci") ;stacksize(2000000) ; 

// for K=2:15, for N=2:5, M(K-l,N-l)=rgindscal2z(N,K) ;end;end;M 

I3=eye(K,K) ; JAC=[] ; IZ= [eye (N-1 ,N-1) , -ones(N-l , 1)] ; 

rankold=0 ; ranknew=l ; 

while ranknew>rankold, 
b=rand(l,N-l) ;bb= [b,-sum(b)] ; c=rand(l ,K) ; 
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JAC= [ JAC ; kron ( IZ , kron (bb , c) ) +kron (bb , kron ( IZ , c ) ) ; kron (bb , kron (bb , 13) ) ] ; 
temp=ranknew;ranknew=rank(JAC) ; 
printf ( 'ranknew=°/oi\n' ,rankiiew) 
raiikold=teiiip; 
end; 

D=ranknew ; 

[lign , col] =size ( JAC) ; Rbar=lign/ (N-l+K) - 1 ; F=Rbar* (N+K-2) -D ; 

function [Rbar,D]=rangjs(N,d) 
// [Rbar,D]=rangjs(N,d) 

// Rbar is the generic rank of sjnninetric tensors of order d and dimension N, 
// and D is the dimension of the whole space of symmetric tensors 
// F=N*Rbar-D is the dimension of the variety of solutions 
// P.Comon, 10 January 1996; modified for scilab 20 may 2006. 
// Call: stacksize (2000000) ; 

// getf ("rangjs . sci") ; f or i=2:6,for j=3:4,M(i-l, j-2)=rangjs(i, j) ;end;end;M 
// NB: preferably limited to d=3 or d=4 and to N<7 for reasons of memory space 
I=eye(N,N) ; JAC=[] ; 
if d==3, 

rankold=0 ; ranknew=l ; 

while ranknew>rankold, 
a=rand (1 , N , "normal " ) ; 

JAC= [JAC ; kr on ( I , kr on ( a , a) ) +kr on ( a , kron ( I , a) ) +kr on (a, kron (a, I))] ; 
temp=ranknew; [Q ,U] =gsm( JAC) ; ranknew=rank (U) ; 
printf ( 'ranknew=°/oi\n' ,ranknew) 
rankold=temp ; 
end; 

D=ranknew ; 
elseif d==4, 
rankold=0 ; ranknew=l ; 
while ranknew>rankold, 
a=rand(l ,N, "normal") ; 

JAC= [JAC ; kr on ( I , kr on ( a , kron ( a , a) ) ) +kr on ( a , kron ( I , kron (a , a) ) ) +kron (a , kron (a , kron ( I , a) ) ) +kron (a , ki 
temp=ranknew; [Q ,U] =gsm( JAC) ; ranknew=rank (U) ; 
printf ( 'ranknew=%i\n' ,ranknew) 
rankold=temp; 
end; 

D=ranknew ; 
else 

Rbar=0;D=0; 
end; 
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[lign, col] =size( JAC) ;Rbar=lign/N-l ; 
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